Supplementary Report for the article Functional diversity can facilitate the collapse of an undesirable ecosystem state

microxanox_version_here <- "0.9.1"
microxanox_version_path <- paste0("microxanox_v", microxanox_version_here)
if (packageVersion("microxanox") < package_version(microxanox_version_here)) {
  stop(paste0("microxanox version needs to be at least ", microxanox_version_here))
}
print(paste0("Version of microxanox: ", microxanox_version_here))

## [1] "Version of microxanox: 0.9.1"

1 Introduction

This document is the Supplementary Report for the report Functional diversity can facilitate the collapse of an undesirable ecosystem state.

2 Growth and inhibition functions

The rate of change in population densities is given by the growth rate - death rate Bush et al 2017. Growth rate is the growth rate in absence of any inhibition multiplied by the amount of inhibition.

Growth rate in the absence of inhibition was modeled with the Monod equation with multiplicative kinetics when two substrates determine growth rate:

\(g_j(X, Y) = g_{max,j} (\frac{X}{K_{j,X} + X})(\frac{Y}{K_{j,Y} + Y})\)

\(g_{max,j}\) is the maximum specific growth rate of strain \(j\). \(X\) is the concentration of substrate 1, and \(Y\) is the concentration of substrate 2. \(K_{j,X}\) and \(K_{j,Y}\) are the half-saturation constants of species \(j\) on substrates \(X\) and \(Y\).

Inhibition of growth by a substrate was modelled with the Haldane equation:

\(h_j(X) = \frac{1}{1+ (X / H_{j,X})}\)

\(H_{j,X}\) can be thought of as a half-inhibition constant: the concentration of \(X\) at which the growth rate of species \(j\) is reduced by 50%.

For further information details of the functions used to model system, please see the original publication by Bush et al. (2018)

3 Creating among strain diversity

Here we provide additional description of the method used to create strain variation within functional groups. The following should be read only after the relevant text in the main report. We generated variation in maximum growth rate \(G_{max}\) and tolerance to reduced sulfur \(h_{SR}\) or oxygen \(h_{O}\) by varying the range of trait values around a reference trait value. To manipulate the amount of variation in a given trait and functional group, we introduced the meta-parameter \(Div\) subscripted by the trait and the functional group. Take for example the cyanobacteria functional group, a reference maximum growth rate parameter \(G_{max,CB}\) and a reference tolerance parameter \(h_{SR,CB}\). The amount of among strain variation in the cyanobacteria is controlled by the two meta-parameters \(Div_{G_{max,CB}}\) and \(Div_{h_{SR,CB}}\). To create a trade-off, the two meta-parameters were always of different sign. The among strain variation was calculated such that the growth rate of strain \(i = 1\) was a factor \(1/({2Div_{G_{max,CB}}})\) of the reference growth rate, and the growth rate of the \(i = n\) strain was a factor \(2Div_{G_{max,CB}}\) of the reference growth rate. The growth rate of the \(i = {2,...n-1}\) other strains was the reference growth rate multiplied by a factor that was equally distributed between the factor of the \(i = 1\) and \(i = n\) strain. This implementation of the variation and trade-off has the assumption of a linear trade-off of log-log transformed among-strain variation. Hence, for example, with the reference growth rate \(G_{max,CB} = 0.05\), among strain variation of \(Div_{G_{max,CB}} = 1\), and nine strains (\(n = 9\)), the growth rates of the nine strains are 0.025, 0.03, 0.035, 0.042, 0.05, 0.059, 0.071, 0.084, 0.1. Put another way, the growth rate parameter of the \(i\)-th of the \(n\) strains in a functional group was calculated as \(G_{max,i} = D_i * G_{max}\) where \(D_i\) is the \(i\)-th element of the set \(D = {2f(x)Div_{G_{max}} |x ∈ N, x ≤ n}\), where \(f(x) = (x − (n + 1)/2)/((n − 1)/2)\).

This procedure ensures that the range of trait values is independent of the number of strains. When there were only two strains, the two strains took the minimum and maximum traits values, as defined by the \(Div\) variable for that trait. When there were three strains, one had the minimum trait value, one the maximum, and one the mean. With nine strains, one had the minimum trait value, one the maximum, one the mean, and three were equally spaced (on log scale) between the minimum and the mean, and the remaining three between the maximum and the mean.

We show three examples of the created variation among strains within functional groups, one for nine strains (Figure 3.1), one for three (Figure 3.2), and one for two (Figure 3.3). In each of these three examples, three graphs show the diversity among strains in the maximum growth rate trait and the tolerance trait, plotted against each other and therefore also showing the trade-off between the two traits.

Notice that in all three cases, the range of trait variation represented by the strains is the same, even though the number of strains differs. This approach allowed us to independently manipulate trait diversity and number of strains. Also note that maximum growth rate of strains is evenly distributed on a log scale, even though it may appear at first site that they are distributed evenly on a linear scale (e.g. Figure 3.2).

All three examples have the following amounts of trait variation: \(Div_{G_{max,CB}}\) = 0.032, \(Div_{h_{SR,CB}}\) = -0.16, \(Div_{G_{max,SB}}\) = 0.095, \(Div_{h_{O,SB}}\) = -1.938, \(Div_{G_{max,PB}}\) = 0.095, \(Div_{h_{O,PB}}\) = -1.938. These are also the maximum amounts of trait variation explored in the simulations.

Figure 3.1: Trait variation with nine strains in each of the functional groups. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%.

Figure 3.2: Trait variation with three strains in each of the functional groups. Note that the three displayed strains are the 1st, 5th, and 9th strains from the nine strain example shown in Figure 3.1. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%.

Figure 3.3: Trait variation with two strains in each of the functional groups. Note that the two displayed strains are the 1st and 9th strains from the nine strain example shown in Figure 3.1. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%.

A consequence of this variation and of the trade-offs is shown in Figure 3.4, which shows the relationship between realised growth rate and the concentration in the environment of the respective inhibiting substance, e.g. reduced sulfur for the cyanobacteria. In these graphs only the least and most tolerant strains are shown. The lines cross due to the trade-off. At high concentrations of the inhibiting substance the more tolerant strain has higher realised growth rate, while at lower concentrations it is the less tolerant strain that has the higher realised growth rate.

Figure 3.4: Response of realised growth rate to variation in concentration of the inhibiting substrate for two strains, namely the most and the least tolerant strain. As always colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains.

4 Method for finding stable states

4.1 General method

We now discuss two approaches for finding how the stable state of the system varies along the gradient of oxygen diffusivity when there exists the possibility of multiple stable states. In the main article we describe and report results of the method that we term the Temporal method.

Here we also describe and give results of what we term the Replication method. In this, we made separate simulations for each of many values of oxygen diffusivity, ran the simulations for long enough to ensure transients had passed, and recorded the state. We ran multiple simulations per value of oxygen diffusivity, and started them with different initial conditions (either low or high total cyanobacterial abundance), so as to check for presence of alternate stable states. This method was used by Bush et al 2017. We term it the replication method since it is achieved by having many independent replicates of the ecosystem, each with a different oxygen diffusivity, and also a separate replicate for each of the different chosen initial conditions.

As is shown in the figures below, the temporal and replication methods result in very similar regions of bistability.

4.1.1 Replication method, no diversity

Figure 4.1: Stable states found with the replication method. The left column of panels refers to the stable state when oxygen diffusivity is at first high (favouring the oxic state). The right column refers to when oxygen diffusivity is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.

4.1.2 Temporal method, no diversity

Figure 4.2: Stable states found with the temporal method. All other aspects are the same as in Figure 4.1.

4.1.3 Both methods overlaid, no diversity

Figure 4.3: Stable states of the two methods overlaid, the temporal method with a dashed line, and the replication method with a solid line.

4.2 Effect of diversity, replication method

4.2.1 Overall effects of diversity

The replication method shows qualitatively the same effects of functional diversity (Figure 4.4) as those revealed by the temporal method. Namely

• Introducing functional diversity in only the cyanobacteria increases the resilience of the oxic state and slightly decreases the resilience of the anoxic state.
• Introducing functional diversity in only the sulfate-reducing bacteria increases the resilience of the anoxic state and reduces the resilience of the oxic state.
• Introducing functional diversity in only the phototrophic sulfur bacteria reduces the resilience of the anoxic state and has no effect on the resilience of the oxic state.
• Introducing functional diversity in more than one functional group reduces, reverses, or does not change the effect of diversity in individual functional groups. For example, introducing functional diversity in all three functional groups reduces the resilience of the oxic state when trait variation is low but increases the resilience of the oxic state when trait variation is high.

Figure 4.4: Replication method results. Left column of panels: the effect of trait diversity (y-axis) on the positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red), which determine the extent of the region of bistability (green shaded region) in the oxygen diffusivity gradient. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. The thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. Right column of panels: The effect size of trait variation on resilience, with positive values on the y-axis indicating that diversity has increased resilience.

4.2.2 An example of stable states

Figure 4.5 shows an example of stable states along the oxygen diffusivity gradient produced by the replication method for finding stable states. The patterns are similar to those produced by the temporal method for finding the stable states (see Figure 4 in the main report), though they differ in one respect: close to the transition back to the state a group dominates, some of the more tolerant strains occur in the stable states found by the temporal method (Figure 4 in the main report), while they do not occur in stable states found by the replication method (Figure 4.5 here). We understand that the appearance of these strains is transient, and the difference in their appearance results from the differences in starting conditions of the two methods. Specifically, in the temporal method, the initial condition at a new level of oxygen diffusivity is the stable state at the previous level of oxygen diffusivity, while in the replication method the starting conditions are identical at each level of oxygen diffusivity. These differences in starting conditions apparently result in longer duration of transients in the temporal method than the replication method.

An interesting feature of these results is that although the position of the oxic to anoxic tipping point is affected by diversity in the sulfate-reducing bacteria (Figure 4.4 and Figure 3 in the main report) there is no evidence of more tolerant strains of sulfate-reducing bacteria being present in the final system state (Figure 4.5). This phenomenon is explained in Section 6.2.

Figure 4.5: Stable states found with the replication method with intermediate levels of diversity in all functional groups. Figure features are otherwise the same as those in Figure 4.1).

4.3 Effect of abundance change event

As mentioned in the main report, to avoid very low strain densities, we added 1 unit of density to every biological state variable every 1,000 hours. We here show that a slightly different method produced the same results. Specifically, we tested the effect of making an abundance floor of 1 unit that was implemented every 1,000 hours. I.e. every 1,000 hours, any abundances that were lower than 1 unit were set to 1 unit. As shown in Figure 4.6, the results with this floor abundance are similar as when 1 was added every 1,000 hours. Note that a coarser gradient of trait variation was used, hence the saw-tooth patterns. A further difference is the large decrease in resilience of the anoxic to oxic transition above moderate levels of diversity in the SB and CB-SB treatments. This pattern can also occur with addition of 1 abundance unit (e.g., see Figure 5.3), and is discussed in the section that concerns this below (section SB diversity effects).

Figure 4.6: Results with an abundance floor. Left column of panels: the effect of trait diversity (y-axis) on the positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red), which determine the extent of the region of bistability (green shaded region) in the oxygen diffusivity gradient. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. The thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. Right column of panels: The effect size of trait variation on resilience, with positive values on the y-axis indicating that diversity has increased resilience, i.e. delayed the shift to the other state.

4.4 Further thoughts

Here we give some further thoughts regarding the comparison of the temporal and replication method for stable state finding.

A particularly notable difference between the two approaches is in the initial condition. In the temporal approach, the initial condition at a new level of oxygen diffusivity is the stable state at the previous level of oxygen diffusivity, while in the approach of Bush et al 2017 the initial conditions are fixed and arbitrary (though this approach is still useful for illustrating the bistability of the system).

The temporal approach is particularly useful when investigating the effects of strain richness. In preliminary simulations using the Bush et al 2017 approach, we found the counterintuitive result that higher strain richness led to lower resilience. We then recalled that we chose to divide the initial total abundance of the cyanobacterial (CB) functional group equally among the strains (termed a substitutive design in biodiversity manipulation experiments). This meant that the most tolerant CB strain in a three-strain simulation had higher initial abundance than the same most tolerant CB strain in a nine-strain simulation. The greater abundance of the most CB tolerant strain made it possible for a three-strain system to develop to be oxic, whereas the nine-strain system developed to be anoxic, despite identical range of trait variation. We could have chosen to do otherwise, for example start with the least tolerant strains of cyanobacteria as more abundant when the system started with high cyanobacterial abundance, and the most tolerant strain as more abundant when the system starts with low cyanobacterial abundance. Using the temporal approach removed the need to specify initial abundances (except the very first one), and therefore results were not so determined by the assumptions used to determine the initial conditions.

5 Effect of simulation duration

5.1 Diversity effects

Recall that the temporal method for finding stable states involves running the simulation for a certain amount of time, termed the “wait time”, at a constant level of oxygen diffusivity. At the end of that wait time the state of the system is recorded, and then the oxygen diffusivity is slightly increased or decreased, and the system then simulated again for the duration of the wait time.

The following figures 5.1-5.5 show the effects of diversity with long to shorter wait times.

Figure 5.1: Effects of diversity with wait time of 1e6 hours (same figure as in main report). (All else as in (Figure 4.4).

Figure 5.2: Effects of diversity with wait time of 1e5 hours. (All else as in (Figure 4.4).

Figure 5.3: Effects of diversity with wait time of 1e4 hours. (All else as in (Figure 4.4).

Figure 5.4: Effects of diversity with wait time of 1e3 hours. (All else as in Figure 4.4).

Figure 5.5: Effects of diversity with wait time of 1e2 hours. (All else as in Figure 4.4).

5.2 Example of transient dynamics at 1e6 hours

The results below show the stable states as a function of oxygen diffusivity for a system with an intermediate level of diversity in all functional groups, with either 1’000’000 hours at each oxygen diffusivity level (Figure 5.6), or 2’000’000 hours (Figure 5.7). This shows that 1’000’000 hours duration of the simulations are sufficient to reach stable state, except very close to state transitions on the trajectory of recovery (i.e.: cyanobacteria: trajectory of increasing oxygen diffusivity; sulfate-reducing bacteria: trajectory of decreasing oxygen diffusivity). Very close to these transitions, even at 2’000’000 hours, the dynamics are still transient, with high tolerance strains at high density, even though they will eventually be outcompeted by low tolerance (higher growth rate) strains.

Figure 5.6: Stable states found with the temporal method with duration at each oxygen diffusivity value of 1e6 hours. The left column of panels refers to the stable state when oxygen diffusivity is at first high (favouring the oxic state). The right column refers to when oxygen diffusivity is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.

Figure 5.7: Stable states found with the temporal method with duration at each oxygen diffusivity value of 2e6 hours. All other aspects are as in Figure 5.6.

6 SB diversity effects

6.1 SB diversity effect on the anoxic to oxic transition

We showed above that wait times of 1e5 (Figure 5.2) and 1e4 (Figure 5.4) produce qualitatively and often quantitatively very similar results to those for wait times of 1e6 (Figure 5.1). There is, however, a notable difference when there is variation in only SB, or in CB-SB with wait times of 1e4 (Figure 5.4). Here, moderate increases in diversity have the same effect as with wait times of 1e6. However, larger increases in diversity suddenly cause a large reduction in resilience of the anoxic-oxic transition. Figure 6.1 gives an example of the stable states when diversity is moderate – increases in oxygen diffusivity cause strain replacement in the SB to the more tolerant strains, which then delay the shift to the oxic state. Suprisingly, when there is just a little more diversity in the SB, then the anoxic to oxic transition happens earlier, i.e. there is a large decrease in resilience (Figure 6.2).

Figure 6.1: Stable state of the system with moderate level of diversity in only SB, and no diversity in the other two groups. All other aspects are as in Figure 5.6.

Figure 6.2: Stable state of the system with higher than moderate level of diversity in only SB, and no diversity in the other two groups. All other aspects are as in Figure 6.1.

We believe that this phenomenon is caused by an effect that we explain with reference to an analogy of using stepping stones to cross a river. Increasing diversity in our simulations is akin to increasing the width of the river, but keeping the number of stepping stones (strains) constant. So long as the river is not too wide, we can still cross it, but once the width of the river is such that the stepping stones are too far apart, we can’t cross the river. And this happens suddenly, the width increases just a tiny bit more, and we suddenly can’t cross the river when we could before.

In the simulated microbial ecosystem, the strains are the stepping stones, and it is harder (takes longer) for strains to replace each other when there is greater distance (in trait space) between them. When they can no longer replace each other in the time available, the system behaves as if only the least tolerant, highest growth rate strain (lightest shade of colour) is present. In the case illustrated above, this leads to much lower resilience of the anoxic system state.

This phenomenon is dependent on having static trait values for the strains. If strain traits were dynamic (e.g. due to evolution), then we expect this result to not occur, and rather for the effects of diversity with two or three strains to be the same as those with six or nine.

6.2 SB diversity effect on the oxic to anoxic transition

When the environment ameliorates for the sulfate-reducing bacteria (i.e. trajectory of decreasing oxygen diffusivity), there is only limited evidence of strain replacement (in final state) along the oxygen diffusivity gradient – either the system is oxic or the least tolerant strain dominates. And yet the switch to anoxic state occurs at higher oxygen diffusivity than when there is no diversity, which suggests there is some role of the more tolerant strains. We here show that the more tolerant strains do play a role, but only a transient one. The following dynamics are for the system starting oxic, and with a value of oxygen diffusivity (log10(oxygen diffusivity) = -4.8). When there is no diversity within functional groups (i.e. two strains that are identical) the system remains oxic (Figure 6.3). When there is diversity (i.e. the two strains differ) then the system switches to the anoxic state (Figure 6.4). In these examples there is only diversity in the two groups of sulfur bacteria. The most tolerant strain does at first grow, but is then outcompeted by less tolerant strains as the concentration of the inhibiting substrate (oxygen) declines.

Figure 6.3: Dynamics of the system with no diversity in any functional group, and log10(oxygen diffusivity) = -4.8. The system remains oxic when starting in the oxic state. Although two strains are shown in the legend, they have identical parameters, and so their dynamics are identical and only one is visible in the dynamics.

Figure 6.4: Dynamics of the system with diversity in the two groups of sulfur bacteria and otherwise exactly the same conditions as in Figure 6.3. The system flips from oxic to anoxic, whereas it did not when there was no diversity.

7 PB diversity effects

Increases in the diversity within the PB functional group leads to lower resilience of the anoxic state (Figure 5.1, Panel PB, blue line). For example, a small amount of variation in only the PB functional group, and no variation in other functional groups, results in transition from anoxic to oxic at approximately \(log_{10}\)(oxygen diffusivity) = -2.76 (Figure 7.1). Whereas a moderate amount of variation in only the PB functional group, and no variation in other functional groups, results in transition from anoxic to oxic at approximately \(log_{10}\)(oxygen diffusivity) = -3.34 (Figure 7.2).

Figure 7.1: Right-hand panels show the transition from anoxic to oxic state when oxygen diffusivity decreases. The results here are for a low amount of variation in only the PB functional group. The transition from anoxic to oxic occurs at approximately \(log_{10}\)(oxygen diffusivity) = -2.76. All other aspects are the same as in Figure 4.1.

Figure 7.2: Right-hand panels show the transition from anoxic to oxic state when oxygen diffusivity decreases. The results here are for a moderate amount of variation in only the PB functional group. The transition from anoxic to oxic occurs at approximately \(log_{10}\)(oxygen diffusivity) = -3.34. All other aspects are the same as in Figure 4.1.

This result is different to the effect of diversity within the CB and SB functional groups, which both cause increases in resilience of the state that they dominate (Figure 5.1). Why does diversity in the PB functional group reduce the resilience of the state it dominates?

The counter-intuitive effect of diversity in the phototrophic sulfur bacteria (PB) possibly results from their positive effects on conditions and processes that favour the oxic state. The phototrophic sulfur bacteria oxidise reduced sulfur to oxidised sulfur. Since reduced sulfur inhibits the cyanobacteria (CB), PB can have a positive effect on CB via their consumption of reduced sulfur. Hence, high growth rates of PB will have a positive effect on CB, and therefore favour the oxic state. We would expect to see that higher diversity of PB is associated with lower concentrations of reduced sulfur and higher concentrations of oxidised sulfur before the transition.

Another pathway by which a higher growth rate of PB can favour the oxic state is by indirectly increasing the concentration of oxygen, thus inhibiting themselves and sulfate-reducing bacteria. This effect occurs because PB consume reduced sulfur, such that less reduced sulfur is available for abiotic oxidation. Consequently, a larger amount of oxygen remains in the environment and inhibits sulfur bacteria. We would expect to see that higher diversity of PB is associated with higher concentrations of oxygen before the transition.

Among these two pathways, effects of phototrophic sulfur bacteria via abiotic oxidation of sulfide are probably of subordinate importance in natural environments because abiotic oxidation rates of sulfide are very small compared to biotic oxidation rates (Luther et al 2011).

There is, however, at least one pathway by which PB can favour the anoxic state: PB produce oxidised sulfur which is consumed by SB, and SB produce reduced sulfur, which inhibits CB.

Looking at the strain replacement that occurs across the oxygen diffusivity gradient, we see that as oxygen diffusivity increases, the more oxygen tolerant strains of PB replace the less tolerant ones because they have higher realised growth rates under the prevailing environmental conditions. Greater PB diversity leads to the availability of strains with higher oxygen tolerance and thus even higher realised growth rates, which via the pathways described above, can favour the oxic state.

Substrate concentrations close to the anoxic-oxic transition in simulations with low and high PB diversity, respectively, corroborate the expected effects (Figure 7.3). As the transition is approached, higher PB diversity leads to higher concentration of oxygen, higher concentration of oxidised sulfur (SO), and lower concentration of reduced sulfur (SR), than a simulation with lower PB diversity. All these differences in substrate concentrations are in the direction of favouring the oxic state, resulting in the earlier transition to the oxic state.

Figure 7.3: Concentrations of three substrates (O: oxygen; SO: oxidised sulfur; SR: reduced sulfur) in a low PB diversity simulation (x-axis) versus substrate concentration in a high PB diversity simulation (y-axis). The figures show only a small range of oxygen diffusivity on the trajectory of increasing oxygen diffusivity. The black line is the 1:1 line.

8 Effect of number of strains

In the main text, we present results of simulations with nine strains per functional group. Here we also explore the effect of number of strains (2, 3, 6, 9) on the position of the two tipping points (Figure 8.1). At low to moderate amounts of trait variation, the number of strains did not affect the position of the tipping points. However, at high amounts of trait variation, resilience was lower when the functional groups had only two or three strains rather than six or nine. This effect was particularly pronounced when there was high variation in only SB or in SB-CB: here, simulations with two or three strains led to dramatically lower resilience of the anoxic state than simulations with six or nine strains. We believe that this effect is due to the greater distance in trait space between strains in simulations with two or three strains than in simulations with more strains (also see Section 5).

Figure 8.1: Effects of number of strains and of diversity among strains.

9 Underlying mechanisms

In the main text we show that trait variation can affect resilience. In this supplementary section, we discuss the mechanism(s) underlying the observed diversity effects. We first explain the two major classes of mechanisms known to drive diversity effects, and then use additional simulations to elucidate the mechanisms operating in our study.

9.1 Complementarity and sampling effect

Research about the relationship between biodiversity and ecosystem functioning has identified two major mechanisms by which diversity affects ecosystem functioning: the selection/sampling effect and the complementarity effect (Loreau & Hector 2001). The sampling effect occurs when more diverse communities have higher functioning because they have a higher probability of containing the best-performing species, i.e. the diversity effect is driven by only one (or few) species (Cardinale et al. 2011). In contrast, the complementarity effect occurs when higher ecosystem functioning in more diverse communities is due to niche differences or positive interactions among species, i.e. the diversity effect is driven by many species (Cardinale et al. 2011). Both effects are considered valid and important mechanisms by which diversity affects ecosystem functioning (Cardinale et al. 2006, 2011, 2012; Fargione et al. 2007).

Similarly, the complementarity effect and the sampling effect are important mechanisms by which diversity affects stability (Loreau & Mazancourt 2013; Loreau et al. 2021) and resilience (Steiner et al. 2006; Hughes & Stachowicz 2011). For example, diversity can stabilize aggregate ecosystem properties (e.g. community biomass) when species differ in their responses to environmental factors and therefore respond asynchronously to environmental fluctuations (temporal complementarity effect; Loreau & Mazancourt 2013). However, when considering ecosystem responses to directional environmental change or pulse disturbances rather than to environmental fluctuations, the sampling effect might be of particular importance (Steiner et al. 2006): more diverse systems might absorb a greater amount of change by having a higher chance of containing species with high environmental tolerance. In this scenario, the positive diversity effect on resilience is driven by one (or few) species. From a conservation perspective, however, one would nevertheless strive to maintain diversity and not only the most tolerant species because of the uncertainties associated with environmental change: we do not know with sufficient certainty which stressors will change in a given ecosystem and which traits will therefore be relevant in the future.

Figure 9.1: A simple phenomenological model of the effect of number of strains in an ecosystem on the resilience of that ecosystem, when the resilience of a ecosystem is linearly proportional to the maximum trait value of the strains present in that ecosystem. Each data point is a different ecosystem, with the strains in each ecosystem drawn at random from a pool of 20 strains. Each of the 20 strains was assigned a trait value drawn at random from a uniform distribution with minimum 0 and maximum 1. The blue line is the fit of a generalised additive model and the grey ribbon shows the 1 standard error confidence range.

9.2 Do the observed diversity effects operate via the sampling effect?

We investigated if the diversity effects observed in our study were in line with the sampling effect, that is, driven entirely by the presence of the most tolerant strain. Specifically, we ran additional simulations where functional groups contained only the strain with greatest tolerance and lowest maximum growth rate. That is, for each of the 20 levels of trait variation, we made simulations where each functional group had only one strain, and its trait values corresponded with those of the most tolerant strain for the given level of trait variation. We also manipulated in which of the functional groups tolerance increased. For example, in a simulation to be compared with variation in only the SB group (and no variation in the CB or PB groups), the SB group contained the strain with highest tolerance, whereas CB and PB each contained the strain with mean tolerance. We then compared the tipping point positions resulting from this tolerance manipulation with the tipping point positions resulting from the diversity manipulation presented in the main text. If the observed diversity effects are caused by the sampling effect, we expect there to be no difference between the tolerance manipulation and the diversity manipulation.

To exclude that differences between the tolerance manipulation and the diversity manipulation were driven by the difference in number of strains, we also performed simulations in which we retained nine strains in each functional group but set the tolerance of all nine strains to the maximum tolerance value. Results were the same as those presented below from simulations with one strain.

For most of the explored combinations and amounts of functional diversity/tolerance the tipping point positions did not differ between the diversity manipulation and the tolerance manipulation (Figure 9.2). However, in some cases there were effects, three of which are described in the follow sections.

Figure 9.2: Effects of strain variation on position of tipping points compared between simulations with all (nine) strains present and simulations where only the most tolerant strain is present.

9.3 High tolerance vs. high diversity in only the SB

When SB contained multiple strains varying in growth rate (and tolerance) there were only two stable states (Figure 9.3), as was previously reported. When, however, only one SB strain with very high tolerance and low maximum growth rate was present, there three stable states: the two previously found, and another where all three functional groups coexisted (Figure 9.4) – this never occurred when there were multiple strains varying in tolerance.

Apparently, at low oxygen diffusivity levels, the absence of SB strains with high growth rate precludes the formation of a stable state that does not include cyanobacteria. It seems that the most tolerant SB strain does not have a high enough growth rate to itself exclude the cyanobacteria. Cyanobacteria are therefore always present, regardless of the oxygen diffusivity. We do not have an explanation of why the cyanobacteria are, however, excluded at intermediate oxygen diffusivities.

Figure 9.3: Stable states found when there is variation among nine strains of SB and the amount of variation is high. All other aspects are the same as in Figure 4.1.

Figure 9.4: Stable states found when SB contained only the most tolerant strain, with trait values corresponding with the diversity level of the previous figure. All other aspects are the same as in Figure 4.1.

In the ecosystem shown in Figure 9.4, where the SB contained one strain with very high tolerance, there was large variation in state variables along the oxygen diffusivity gradient at values less than about -4.5. To examine what might be causing this, we simulated the ecosystem at a constant oxygen diffusivity of \(log_{10}\) oxygen diffusivity of -6 (Figure 9.5). We see that the dynamics are qualitatively different from other cases in which there is a stable state reached in such a simulation (Figure 6.4). It appears that the high tolerance SB strain is able to initiate a transition to an anoxic state, but is not able to create a stable anoxic state. This dynamical situation was not present in the ecosystem in which high tolerance and low tolerance strains were present.

Figure 9.5: Temporal dynamics of an ecosystem with a single SB strain with high tolerance, the same one as in Figure 9.4 simulated with constant \(log_{10}\) oxygen diffusivity of -6.

9.4 High tolerance vs. high diversity in CB-SB

In the SB-SB case, presence of only the high tolerance strain could cause the system to not display any alternate stable states (Figure 9.6).

Figure 9.6: Stable states found when CB and SB each contain only the single most tolerant strain with trait values corresponding with the diversity level of the previous figure. All other aspects are the same as in Figure 4.1.

9.5 High tolerance vs. high diversity in CB-SB-PB

When in each of the three functional groups only the most tolerant strain is present, the resilience of the oxic state is greater than when there is variation in strain tolerance within groups. It is not entirely clear why this effect only occurs only in this combination.

Figure 9.7: Stable states found with nine strains of CB-SB-PB and high diversity. All other aspects are the same as in Figure 4.1.

Figure 9.8: Stable states found with the single most tolerant CB strain, SB strain, and PB strain, with trait values corresponding with the diversity level in the previous figure. All other aspects are the same as in Figure 4.1.

9.6 Why only at high diversity/tolerance?

Why do these effects only occur at high diversity/tolerance? The likely reason is that at high diversities the most tolerant strain has particularly low growth rate, and so it may have the tolerance required to increase and persist, but not the growth rate required to dominate.

10 Finding the separatrix

In the original publication (Bush et al 2017) Figure 3a shows a separatrix in the cyanobacteria abundance - oxygen diffusivity dimensions of the system. If the system has initial cyanobacterial abundance greater than the value at the separatrix then the cyanobacterial abundance increases, and the system goes to the oxic stable state. Conversely, if the system has initial cyanobacterial abundance lower than the value at the separatrix then the cyanobacterial abundance decreases, and the system goes to the anoxic stable state.

It is critical to note, however, that the position of the separatrix depends on both the parameters of the system, and the initial conditions of variables other than cyanobacterial abundance. If the system were started with higher initial concentration, then the position of the separatrix would change (as would the values of oxygen diffusivity at which the system tipped from anoxic to oxic, and from oxic to anoxic).

This issue of dependence of the position of the separatrix on the value of initial condition of multiple state variables has particular significance when there are multiple strains per functional group. Then one cannot just vary initial total cyanobacterial abundance (as was done to create the separatrix on Figure 3a in the original publication) but one must also define how the initial abundance of each of the strains varies. One might choose to make each strain have equal initial abundance, and to vary the total (as we did in our implementation of the replication method in the main article).

We could, however, and arguably should, give the most tolerant strain of cyanobacteria higher initial relative abundance when total initial cyanobacterial abundance was low, and the least tolerant strain higher initial relative abundance when total cyanobacterial abundance was high. And we would need to specify the total and relative abundance of all cyanobacteria (and sulfur bactria) strains. Thus, with a system as complex as the 9 strain system, which has 31 state variables the separatrix is very high dimensional, and evaluating (and plotting) it in one dimension requires many assumptions. We chose to not attempt to find a meaningful set of assumptions in this study, though not because we think it would be unprofitable to do so, but rather because we did not want to risk distracting attention from findings we consider to be significant enough in their own right.

11 (Non-)additivity

We investigated if the effects of trait variation of individual functional groups were additive. We compared the tipping point positions as determined by our simulations with the expected position of the respective tipping point if effects of diversity in different functional groups were additive. The expected tipping point was calculated by adding the effect sizes (i.e. the difference in tipping point position between the no-diversity scenario and a given amount of trait variation) on a linear scale.

Figure 11.1: Assessing the additivity or non-additivity of the effects of trait variation. Solid black lines depict the position of the two tipping points at different levels of trait variation, as simulated with the model. The colored, dashed lines show the expected position of tipping points if effects of diversity in different functional groups were additive. The expected tipping point was calculated by adding the effect sizes of trait variation of different functional groups on a linear scale. For example, in (d) the effect size of variation in only CB was added to the effect size of variation in only SB. In (g), the effect size of variation in only CB was added to the effect size of simultaneous variation in both SB and PB. (CB: cyanobacteria, SB: sulfate-reducing bacteria, PB: phototrophic sulfur bacteria). The dotted vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations.

When two functional groups varied in traits, the effects of trait variation were additive when one of the groups had no diversity effect on a tipping point but non-additive when both groups had diversity effects on a tipping point (Fig. 11.1 a-f). High trait variation in cyanobacteria and sulfate-reducing bacteria caused the oxic-anoxic tipping point to occur at lower oxygen diffusivity than expected assuming additivity of the individual effects of diversity in these two groups (Fig. 11.1d). When the two groups of sulfur bacteria varied in traits, the anoxic-oxic tipping point occurred at higher oxygen diffusivity than expected by additivity (Fig. 11.1f).

When all three groups varied in traits, only the diversity effect in phototrophic sulfur bacteria on the oxic-anoxic tipping point was additive; all other effects were non-additive, especially at higher trait variation (Fig. 11.1g-i).

12 Comparing the effects of different treatments

Figure 12.1: Another method for comparing the effects of different diversity treatments. Some lines (treatment combinations) are over-plotted and therefore show identical effects of some combinations.

13 CB diversity effects

Increased CB diversity causes a slight decrease in the resilience of the anoxic state. The following two figures show the state responses to oxygen diffusivity with low Figure 13.1 and high Figure 13.2 levels of CB diversity.

The difference in the position of the anoxic to oxic tipping point is small. It is, in fact the smallest it could be, as it is the smallest step size in oxygen diffusivity used in the simulations. Furthermore, we think the presence of higher strain diversity can increase the rate of response of the system (since there is a strain with higher growth rate than in low diversity levels) and that this difference in response rates is likely responsible for the very small difference in resilience.

Figure 13.1: Stable state of the system with low level of diversity in only CB, and no diversity in the other two groups. All other aspects are as in Figure 5.6.

Figure 13.2: Stable state of the system with high level of diversity in only CB, and no diversity in the other two groups. All other aspects are as in Figure 5.6.

14 References

Bush, T., Diao, M., Allen, R.J., Sinnige, R., Muyzer, G. & Huisman, J. (2017). Oxic-anoxic regime shifts mediated by feedbacks between biogeochemical processes and microbial community dynamics. Nat. Commun., 8, 789.

Luther, G.W., Findlay, A., MacDonald, D., Owings, S., Hanson, T., Beinart, R., et al. (2011). Thermodynamics and kinetics of sulfide oxidation by oxygen: a look at inorganically controlled reactions and biologically mediated processes in the environment. Front. Microbiol., 2, 62.